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A general analytical formula for recurrence relations of multisite interaction Ising mod- 
els in an external magnetic field on the Cayley-type lattices is derived. Using the theory 
q ■ of complex analytical dynamics on the Riemann sphere, a numerical algorithm to ob- 

tain Yang-Lee and Fisher zeros of the models is developed. It is shown that the sets 
of Yang-Lee and Fisher zeros are almost always fractals, that could be associated with 
Mandelbrot-like sets on the complex magnetic field and temperature planes respectively. 

O ! 1 Introduction 

In 1952, Lee and Yang in their famous papers JTJ first proposed a new and effective method 
for investigation of phase transitions. They studied the distribution of zeros of the partition 
function considered as a function of a complex magnetic field (fugacity), and proved the circle 
theorem, which states that the zeros of the partition function of Ising ferromagnet lie on the 
unit circle in the complex fugacity plane ( Yang-Lee zeros) . After these pioneer works of Lee and 
Yang, Fisher [Q, in 1964, initiated the study of partition function zeros in the complex tem- 
perature plane (Fisher zeros). These methods are then extended to other type of interactions 
and were widely applied || . Recently, Binek j| showed that the density of Yang-Lee zeros can 
be determined experimentally from the field dependence of the isothermal magnetization data. 
rN [ The fractal structure of Fisher zeros in q-state Potts model on the diamond lattice was obtained 
by Derrida, Seze and Itzykson ||. They showed that the Fisher zeros in hierarchical lattice 
models are just the Julia set corresponding to the renormalization transformation. Bosco and 
Goulard Rosa |J investigated Yang-Lee zeros of the ferromagnetic Ising Model on the Cayley 
tree and associated these with the Julia set of the Cayley tree recursion map. They studied the 
case when the complex magnetic field was applied only to spins on lattice surface and found 
that Yang-Lee zeros are distributed along the unit circle. 

The lattice models with multisite interaction Ising and Heisenberg models have been used 
for modelling the physical systems such as the binary alloys ||, classical fluids ||, liquid bi- 
layers [[LI]], solid 3 He [fTU| , rare gases Jl2| and anisotropic magnets (CeBi, EuSe, etc.). These 



systems have rather complicated phase diagrams and unusual properties. Indeed, the Monte 
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Carlo results give evidence of the presence of phase transitions at nonzero values of the magnetic 
field for a number of ferromagnetic multisite interaction models ||13|| . There are few analytical 
results for these systems due to their greater complexity. Here we proposed a numerical algo- 
rithm for investigation of complex zeros of partition functions of models on the Cayley-type 
lattices. 

This Letter is organized as follows. In Section 2 we give the model and derive the generalized 
recurrence equation. In Section 3 it is shown that the zeros of partition function can be 
associated with the set of external parameters (kT, magnetic field) at which the recurrence 
function has neutral periodic cycles. In Section 4 we present some definitions and theorems 
from complex analytic dynamics and give the numerical algorithm for obtaining Yang-Lee and 
Fisher zeros of the models. In the last section we give the explanation of the results and make 
conclusions. 



2 The Model 

We consider a multisite interaction Ising model on the Cayley-type lattice. The Cayley 
tree and Husimi lattice are the well known representatives of this class of recurrence lattices. 
In general, the Cayley-type lattice is constructed of p-polygons. It is characterized by p, the 
number of edges (the number of sites) of the polygon (p = 2 - usual Cayley tree, p = 3 - Husimi 
lattice) and by q, the number of p-polygons that go out from each site (Fig.l). 




Fig. 1. The recursive structure of 4-polygon Cayley-type lattice with q = 3. The 
numbers 0,1,2 stand for shells. 

One of the pecular properties of the Cayley-type lattice is that in the thermodynamical 
limit (the number of shells are tended to oo) the number of surface sites becomes proportional 
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to the number of inner sites of the lattice (for more detailes see [14 , [T 



Therefore, the models 
In the infinite 
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defined on the Cayley-type lattice exhibit quite unusual behavior 
dimensional Euclidean space the Cayley-type lattices are usually viewed as ramified trees with 
constant vertex connectivity. However, these can be embedded in the two-dimensional space of 
constant negative curvature (the hyperbolic plane) with fixed bond angles and length M. 

The Hamiltonian of the multisite interaction Ising model on a p-polygon Cayley-type lattice 
has the form 

(polygons) i 

where Si takes on values ±1, the first sum goes over all the p-polygons of the lattice and fl is 
the product of all spins on a p-polygon; the second sum goes over all the sites on the lattice. 

The Cayley-type lattice has the advantage that for models defined on it one can derive the 
exact recurrence relation and use the theory of dynamical systems for investigation thermody- 
namical properties of the models fl4]]. Let us denote the statistical weight of the n-th generation 
branch with the basic spin at the state S as g n {S). By cutting of the n-th generation branch 
at the basic p-polygon one obtains (q — — 1) interacting (n — l)-th generation branches 
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w(S , Si,..., Sp-^gKSi) ■ ■ ■ gl(S p -i), 



(2) 



where S , Si, ... , S p -i are the spins on the basic p-polygon, 7 = q — 1 and w(S , Si, ... , S p _i) 
is the statistical weight of the basic p-polygon 



w(S , 5p_i) = exp 



p-i 1 p-i 

jp IT h — Sk 

k=0 1 k=0 



(3) 



where J p = J p /kT, h = h'/kT. Cutting apart the lattice at the central site one can obtain for 
the partition function Z = J2 e~@ H the following expression 



Z n = J2gl(S ). 

So 



(4) 



Using Eq.(|) and introducing an auxiliary quantity x n = g n (+)/gn(~), we can get the recur- 
rence relation for x n , and express the thermodynamic quantities such as the magnetization, the 
specific heat, the free energy etc., in terms of x n [0 . 
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J2k=0 O-k+lCp^X^ 1 



2^k=0 U-k^p-lX 



where 




exp 
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(5) 



(6) 



and C k = p\/(k\(p — k)\) is the simple combinations of p elements taken A; at a time. Using the 
binomial theorem, one can obtain after some calculations the following general formula for the 
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recurrence relation of the multisite interaction Ising model on p-polygon Cayley-type lattice in 
an external magnetic field x n+ \ = f(x n ), where 

_ (e 2h x q - 1 + l)^ 1 + tanh J v {e 2h x^~ x - Vf' 1 

~ \ e 2h x q-l + _ tanh J p l e 2h x q-1 _ ^p-V ^ 

2h 

Here we used the substitution Xk — > e « Xk- 

In the following we shall use some definitions from the theory of dynamical systems. For 
every xq on the Riemann sphere (xq G C) x n+ i = f(x n ), n = 0, 1, 2, . . . , generates the forward 
orbit of Xq, which is denoted by Or + (x )- If x n = x* for some n in Or + (x*) we say that x* is 
a periodic point, and Or + (x*) is called a periodic orbit or periodic cycle. If n is the smallest 
integer with that property, then n is the period of the cycle. Usually, if n — 1, x* is called a 
fixed point. Obviously, x* is a fixed point of f n , if x* is a periodic point of period n. (One 
should not confuse the iterations of / with the powers of /, i.e. f n = f o . . . o / is different from 
[f(x)] n ). To characterize the stability of a periodic point x* of period n we have to compute the 
derivatives. The complex number A = (f n )'(x*) = -^j is called the eigenvalue of x*. Using 
the chain rule of differentiation we see that this number is the same for each point in a cycle. 
A periodic point x* is called (1) attracting (stable) if |A| < 1, (superattracting if A = 0), (2) 
repelling if |A| > 1, (3) rationally neutral (indifferent) if |A| = 1 and A n = 1 for some integer n, 
(4) irrationally neutral if |A| = 1, but A n is never 1. 



3 Zeros of the Partition Function and Neutral Periodic 
Cycles 



Let us now consider the partition function of the model 

Z n ~gl{-){e 2h xl + l)- > 



9 q nH(e 2h [f n (x )r + l), 



(8) 



where g n {~) is an analytic function and g n {—) 7^ 0, oo for any kT, h and n; Xq is an initial 
point of iterations (x = 1 corresponds to the free boundary condition). From Eq.(|8|) one can 
see that Z n is a rational function and the free energy of the system has the form 



-kT lim In Z„ 



-kT lim In 



9 q n( 



Qn(z,fx) 



(9) 



where P n (z,fi) and Q n (z,fi) are polynomials in z = e 2Jp and fi = e 2h , where all coefficients are 
positive (see Eq.@). Since g n (—) is an analytic function, the free energy T has a singularity if 
and only if P n {z,\i) = or Q n (z,fi) = 0. It is obvious from Eq.flD and Eq.@ that Q n (z,fi) is 
the denominator of f n {xo) function. Since 00 cannot be a periodical point for f(x) function, the 
condition Q n (z, /i) = will not produce singularity of the free energy T in the thermodynamical 
limit. Hence, the singularities of T (phase transition points) correspond to zeros of the partition 
function Z n (P n (z,fi) = 0), i.e. Yang-Lee or Fisher zeros |1|, 

Let us now consider the problem of phase transitions on the Cayley-type lattices in terms 
of recurrence relations and dynamical systems theory. At high temperatures T > T* (T* is the 
paramagnetic phase transition point) the recurrence relation has only one stable fixed point 
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x*, corresponding to the paramagnetic phase. When the temperature is lowered, two different 
scenarios occur depending on the type of interactions. For ferromagnetic interactions (J p > 0) 
at T < T* the "paramagnetic" fixed point either becomes unstable and there arise two different 
attracting fixed points, each corresponding to one of the two ferromagnetic phases with opposite 
directions of magnetization or the "paramagnetic" fixed point remains stable and there 
arise an additional attracting fixed point (See Fig. 2 (d)). For antiferromagnetic interactions 
(J p < 0) the "paramagnetic" fixed point at T < T* becomes unstable and there occur a period 



doubling bifurcation cascade and even chaos [19]. One says that for given external parameters 



(kT, h) the system is in the stable equilibrium state when the iterations of f(x) approach 
a stable (attracting) fc-cycle, and the system undergoes phase transition when, by changing 
the external parameters, the iterations of f(x) pass from one (nonstable) A;-cycle to the other 
(stable) fc'-cycle (phase transitions between the modulated phases p0[). The values of external 
parameters at which the phase transition occurs may be obtained from the following conditions 

/*(*)=*, \(f)'(x)\ = l. 

We see that the phase transition points could be associated with neutral periodic cycles of the 
mapping. We call the set of parameter values, at which the rational function / from Eq. (0) has 
neutral periodic cycles, the Mandelbrot-like set of /. Thus, to investigate the Yang-Lee and 
Fisher zeros one can study the Mandelbrot-like sets of /. 



4 Mathematical Background and the Numerical Algo- 
rithm 

In this section we give the required definitions and theorems from complex analytic dynamics 
on the Riemann sphere pi], and develop a numerical algorithm for obtaining the Yang-Lee 
and Fisher zeros of the model. 

One can consider any rational function f(x) = P(x)/Q(x), where P(x) and Q(x) are poly- 
nomials in x as an holomorphic map on a Riemann sphere. The poles of the rational function 
are simply the points of Riemann sphere C = C\J {oo}, that are mapped to oo (the upper pole 
of the Riemann sphere). The degree of f(x) d = deg(f) is defined as 

d = max{deg(P), deg(Q)}. 

Also, the degree of / is the number (counted with multiplicity) of inverse images of any point 
of C. v is a critical value of / if the equation f(c) = v has a solution the multiplicity of 
which is greater than one. Such a solution c is called a critical point. In local coordinates this 
is equivalent to the condition f'(c) = (at least when c ^ oo). To discuss the behavior of 
f(x) near oo one usually invokes another transformation, the reflection r(x) = 1/x at the unit 
circle, which exchanges with oo; i.e. to treat the behavior of f(x) near oo, one is to consider 
<p(x) = r(f(r(x))) = roRor(x) in the vicinity of 0. Below we give some theorems from complex 
analytic dynamics without proofs. The reader should consult Carleson and Gamelin |H[ for 
proofs. 

Theorem 1. The number of critical points of / is at most 2d — 2. 

If x* is an attracting fixed point for /, we define the basin of attraction of x*, denoted 
by A(x*), to consist of all x such that f n (x*) is defined for all n > 1 and f n (x) — > x*. The 
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connected component of A(x*) containing x* is called the immediate basin of attraction of x* 
and is denoted by A*(x*). If 7 is an attracting cycle of period n, then each of the fixed points 
of f n (x*) has its basin and A (7) is simply the union of these basins. 

Theorem 2. If x* is an attracting periodic point, then the immediate basin of attraction 
A*(x*) contains at least one critical point. 

Theorem 3. If x* is a rationally neutral periodic point, then each immediate basin of 
attraction associated with the cycle of x* contains a critical point. 

Theorem 4. The total number of attracting and neutral cycles is at most 2d — 2. 

The numerical algorithm is based on Theorems 1-4 and the well known fact that the con- 
vergence of iterations to neutral periodic cycle is very weak; i.e. one has to make a number 
of iterations (infinite number of iterations) in order to approach the neutral periodic cycle. 
Of course, the initial point does not belong to the periodic cycle or the Julia set of mapping 
(for the definition of the Julia set see Ref. The algorithm is as follows: One finds all 

the critical points of the mapping and investigates the convergence of all the orbits of critical 
points (critical orbits) to any attracting periodic cycle. If all critical orbits converge, for ex- 
ample, after n iterations, one says that the system is in stable equilibrium state, otherwise the 
system undergoes phase transition. Of course, the last statement is not rigorous since a weak 
convergence to an attracting periodic cycle is possible. Depending on the choice of n and e 
- the accuracy of approaching the attracting cycle, the phase transition pictures may change. 
Our complutergraphical experiments showed that n = 700 and e = 10 -4 are optimal values for 
n and e, and the pictures generated by this algorithm will not qualitatively change when n is 
increased and/or e decreased from their optimal values. In our algorithm we supposed that the 
cases where neutral periodic cycles are irrational are very rare and negligible. 

5 Results and Pictures 

One can easily find all the critical points of mapping / from Eq.(^) as is described in the 
previous section. The critical points are as follows: x = with multiplicity q — 2, x = 00 
with multiplicity q — 2 and 2(q — l)(p — 1) solutions of equations (e^x 9-1 + = and 

(e 2h x q ~ l — l) p_1 = 0. The degree of our mapping / is d = (p — l)(q — 1) and, according to 
Theorem 1 of Section 4, these all are critical points of /. One can easily find that the most of the 
critical orbits after the first iteration are intersected and can consider only the orbits starting at 
points 0, z (00 — > z), 1, — 1 for p > 2 and 1/z, z for p = 2. Also, the orbits of and 00 intersect 
when p is odd and the orbits of 1 and -1 intersect when q is odd. Bellow we give some pictures 
generated by our algorithm. Fig. 2 shows computergraphical study for Fisher zeros of the Ising 
model on the Cayley tree (p = 2, q = 3) for J' = 1 and different values of magnetic field hi . We 
have experimental evidence that all critical orbits converge in white regions. Physical phase 
transitions take place at the points where black regions intersect the real axis. We see that 
for h! = (Fig. 2 (a)) the second order phase transition occurs for z = 3 in accordance with 
Baxter's result pl| and for h! = 1 (Fig. 2 (c)) there is no phase transition for kT > 0. For 
hi = 0.5 (Fig. 2 (b)) the first order phase transition occurs in "cat face" -like region. Figure 
2 (d) shows the neutral fixed point of f(x) function corresponding to phase transition in the 
" cat face" -like region of Fig. 2 (b). A care should be taken in calculations for small kT < 0.07 
because one of the attracting fixed points tends to infinity when the temperature is lowered. 
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Since the free energy of the Ising model on the Cayley tree is an even function of h', Figs. 2 
(a)-(c) are invariant under the transformation hi — > —hi and show phase transitions for both 
ferromagnetic and antiferromagnetic interactions (kT — > —kT is equivalent to J' — > —J' and 
hi — > — /&'). In Fig. 3 we present Yang-Lee zeros of the ferromagnetic Ising model on the Cayley 
tree. Figs. 3 (a)-(b) show that the phase transition exists at positive temperatures for hi = 0.5 
and hi = —0.5 in accordance with Fig. 2(b). The set of Yang-Lee zeros at the critical point 
z = 3 is shown on Fig. 3(c). One can see that in this case a real phase transition occurs at 
zero magnetic field only as in usual ferromagnetic systems. It is interesting to note that the set 
of Yang-Lee zeros of ferromagnetic Ising model on the Cayley tree resembles the boundary of 
the Mandelbrot set of the quadratic mapping z — > z 2 + c. In this case the Yang-Lee zeros are 
not located on the unit circle because in the thermodynamical limit the contribution of surface 
and inner spins to the partition function of the model on Cayley tree are of the same order. In 
Fig. 3 (d) we give the improved set of Fisher zeros of Ising model on the Husimi lattice (p = 3, 



5 = 4) previously obtained in Ref. P3 



In conclusion, we should like to note that the numerical algorithm presented in this Letter 
may be applied for obtaining Yang-Lee and Fisher zeros of any model on the Cayley-type lattice, 
for which one-dimensional recurrence relation can be derived. It is remarkable that there appear 
the well-known Mandelbrot set figures. This phenomenon is known as the universality of the 
Mandelbrot set p4|. 
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Figure Caption 



Figure 1. The recursive structure of 4-polygon Cayley-type lattice with q = 3. The numbers 
0, 1, 2 stand for shells. 

Figure 2. The Fisher zeros of Ising model on the Cayley tree, (a) J' = 1, h' = 0; (b) J' = 1, 
h! = 0.5; (c) J' = 1, h' = 1; (d) The plot of f(x) function from Eq. (|7|) corresponds to the 
phase transition point kT = 0.7256, J' — 1, h! — 0.5 of (b). 

Figure 3. (a)-(c) The Yang-Lee zeros of Ising models on the Cayley tree; (a) p = 2, q — 3, 
z = 15, 742; (b) a close up from left; it is noteworthy that the set in question may be obtained 
by mapping r(/i) = from the Mandelbrot set boundary of (a); (c) p = 2, q = 3, z = 3; (d) 
The Fisher zeros of multisite interaction Ising model on Husimi lattice J' — — 1, h' — 3, p — 3, 
q = 4. 
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